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We utilize tlie novel non-Markovian quantum jump (NMQJ) approach to stochastically simulate exciton dy- 
namics derived from a time-convolutionless master equation. For relevant parameters and time scales, the 
time-dependent, oscillatory decoherence rates can have negative regions, a signature of non-Markovian behav- 
ior and of the revival of coherences. This can lead to non-Markovian population beatings for a dimer system 
at room temperature. We show that strong exciton-phonon coupling to low frequency modes can considerably 
modify transport properties. We observe increased exciton transport, which can be seen as an extension of recent 
environment-assisted quantum transport (ENAQT) concepts to the non-Markovian regime. Within the NMQJ 
method, the Fenna-Matthew-Olson protein is investigated as a prototype for larger photosynthetic complexes. 
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I. INTRODUCTION 

The initial step in photosynthesis is the excitonic transport 
of the energy captured from photons to a reaction center HI. 
In this process, highly efficient transport occurs between in- 
teracting chlorophyll molecules embedded in a solvent and/or 
a protein environment Q . The exciton transfer dynamics has 
been studied utilizing Forster theory in the limit of weak inter- 
molecular coupling |3| or Redfield master equations in the 
limit of weak exciton-phonon coupling |4|. The latter ap- 
proach describes the transport as dissipative dynamics for the 
reduced excitonic density matrix. Master equations are devel- 
oped starting from projector operator techniques that separate 
relevant (system) from less relevant (phonon) degrees of free- 
dom. Formally exact dynamics for the system is described by 
the Nakajima-Zwanzig equation with time-convolution 1 5 , 6 1 
and the time-convolutionless (TCL) equation [7.,8, 9. 10 .llj . 
The first is equivalent to the chronological ordering prescrip- 
tion (COP) while the second corresponds to a partial ordering 
prescription (POP) of the time-ordering in a system-bath cu- 
mulant expansion ||2l[l2l[I3. The convolution kernel can be 
transformed into the TCL kernel by including the appropriate 
backward propagation for the density matrix ifTOl fT4l . The 
time-dependent Redfield equation is derived from a second- 
order approximation in the system-bath interaction Hamilto- 
nian |14|. Further imposing the Markov approximation leads 
to the standard time-independent Redfield equation. The dy- 
namics of the populations of the density matrix and that of 
the coherences is separated when the secular approximation is 
employed, in which the master equation can be cast into Lind- 
blad form. Recently, Palmieri et al. developed a prescription 
for reintroducing the coupling of populations and coherences 
with suitably defined Lindblad operators Ill5l . 

Non-Markovian (NM) effects can be taken into account 
by the time-convolutionless approach. Other frequently used 
methods explicitly include strongly coupled modes or envi- 
ronmental two-level systems into the system dynamics lfT6l 
[TtI . The Hilbert space size and thus numerical effort increases 
concomitantly. Kubo, Tanimura et al. developed a hierarchi- 
cal treatment where auxiliary systems describe higher order 
system-bath interactions ifTsl [191 l20l . Xu et al. introduced an 
elegant filtering method for this approach ETl 12211 . While the 



hierarchical treatment is formally exact for Gaussian fluctua- 
tions, the infinite set of equations is truncated for numerical 
propagation. The method was recently applied to investigate 
long-lived coherence in the Fenna-Matthews-Olson (FMO) 
protein complex |23 , 24 1. 

A Markovian master equation in Lindblad form can be sim- 
ulated by means of the Monte-Carlo wavefunction method 
(MCWF) ||25l . This numerical technique relies on the prop- 
erty that density matrix evolution is equivalent to an averaging 
of wavefunction trajectories, each of which is interrupted by 
stochastic, discontinuous quantum jumps. In this work, we 
employ the non-Markovian quantum jump approach (NMQJ), 
recently developed by Piilo et al. [26. 27 1. This method is 
a generalization of the MCWF to the case of non-Markovian 
dynamics derived from a TCL approach. The TCL approach 
can lead to time-dependent, oscillatory decoherence rates that 
have negative regions, a signature of NM behavior. These neg- 
ative rates are shown to lead to a reversal of decoherence by 
defining appropriate quantum jumps and jump probabilities. 

Compared to the explicit numerical integration of the mas- 
ter equation, the NMQJ approach has interesting features in 
the context of exciton transfer in chromophoric networks. 
First, the NMQJ approach is, similar to the MCWF, based 
on the propagation of wavefunctions and thus scales consid- 
erably better with the system size than approaches in Liou- 
ville space. It is therefore especially suitable for simulat- 
ing large chromophoric networks of photosynthetic antenna 
systems. Second, positivity violation ifTOl can be efficiently 
detected during the simulation by a simple criterion for the 
negative jump probabilities 1281. Third, the trajectory pic- 
ture allows for new insights into exciton dynamics. Quantum 
jumps related to negative transition rates can restore coher- 
ence and thus can provide an additional theoretical perspec- 
tive on long-lived quantum coherence found in photosynthetic 
systems such as the Fenna-Matthews-Olson complex |23 1 and 
the reaction center of purple bacteria j29|. Here, we apply the 
NMQJ approach to dimer systems and the FMO complex. We 
observe population beatings arising from recurrence effects of 
the NM environment. We also find that in the non-Markovian 
regime transport can be enhanced compared to purely Marko- 
vian dynamics and thereby provide an extension to the recent 
ENAQT concept Il30ll3ni32l . These effects are pronounced in 
situations when the main phonon-mode frequencies are much 



2 



smaller than (i.e. "off-resonant" to) a particular system transi- 
tion frequency. 

In Sec. 2, we develop a time-convolutionless master equa- 
tion for the dynamics of a single excitation, leading to time- 
dependent rates. In Sec. 3, we discuss the spectral density and 
time-dependent rates in more detail and explain the physical 
situation where NM effects are considerable. In Sec. 4, we in- 
troduce the NMQJ method in the context of excitonic energy 
transfer. Finally, in Sees. 5-7, we analyze dimer systems and 
the Fenna-Matthews-Olson complex. 



II. MASTER EQUATION 

The transport dynamics of a single excitation is described 
by a master equation for the density matrix that includes co- 
herent evolution, relaxation, and dephasing. In this work, we 
are mainly interested in the effect of slow phonon fluctuations 
and the memory of the bath on the excitonic energy transport 
dynamics. We utilize a time-convolutionless master equation 
to second order in the exciton-phonon coupling. We employ 
the secular approximation and focus on non-Markovian de- 
coherence rates. The removal of the secular approximation 
requires modifications to the NMQJ approach and is left for 
future work. The validity and limitation of the Redfield ap- 
proach with respect to parameters such as environmental cou- 
pling and temperature and the with repect to the neglect of 
fundamental processes such as the molecular reorganization 
has been discussed in detail in |4, 33 34 1. The complete 
Hamiltonian for an interacting A^-chromophoric system in the 
single exciton manifold and including the phonon part is given 
by iJ = Hs + Hs B + Hb ■ The system part is in tight-binding 
form: 

N N 

Hs=Y^ e™|m)(TO| + Knn(|m)(n| + |n)(m|), (1) 

m—l n<m 

where the Hilbert space basis states \m) denote the presence 
of an excitation at the mth chromophore and e,„ are relative 
site energies with respect to the chromophore with the low- 
est absorption energy. The Vmn are the interchromophoric 
couplings. The exciton basis \M) = X^m ''ni(-^'^)l'™) '■^^ 
eigenbasis of the Hamiltonian Hs\M) = Em\M). The 
exciton-phonon Hamiltonian is dominated by site energy fluc- 
tuations: 

HsB^^A^^B^, (2) 

m 

with the system part Am = \m){m\ and the bath part B„i = 
( Xj /itJiAi(6i + 6|) ) . Each site is separately interacting 

with a set of phonon modes indicated by the subindex m of 
the bath part. The dimensionless coefficients A; describe the 
coupling strength to each phonon mode. The phonon Hamil- 
tonian is Hb = J2i m [fi'^ibibi) , where the sum is over 
all phonon modes (at each site) described by the bosonic op- 
erators bi and frequencies uji. For this work, we assume that 



the chromophores are coupled independently to their respec- 
tive baths. Recent studies include spatial correlations into the 
exciton dynamics ||35l [36l l37ll . 

The second order TCL master equation for the reduced sys- 
tem density matrix in the interaction picture is given by Iil4il : 

JtPsit) = -^J dtitrB{[HsB{t),[HsB{h),p{t)(g>pB]]}. 

(3) 

The characteristic double commutator arises from a second- 
order perturbation treatment of the system-bath Hamiltonian. 
Note that it is assumed that the effect of the system on the 
bath is small such that the total system approximately fac- 
torizes for all times. Multiphonon processes, arising from 
higher order commutators, are not taken into account. An ad- 
ditional approximation in standard Redfield theory is that the 
phonon bath is always in equilibrium, which neglects molecu- 
lar reorganization effects. Next, we introduce the operators 
AM = E.„-.„=a.<„(M)c™(iV)|A/)(iV|, which de- 
scribe the effect of the system-bath Hamiltonian in the eigen- 
basis of the system Hamiltonian, i.e. Am = Am{^^), 
where the sum is over all transitions in the single exciton man- 
ifold fnil . This leads to the master equation of the form: 

|psW = -^^5][A„M,[A„(c.'),pW]] (4) 
Jo 

with the symmetrized correlation function, 

Sm{t) - \tr B{{Bmit),BmmPB}, (5) 

and the associated response function, 

Xmit) = ~i trB{[Bm{t),BmmPB}- (6) 

The quantities in Eq. (j5]l and Eq. (j6]l are related to the real 
and imaginary parts of the bath correlator, i.e. Cm{t) = 
trB{Bmit)Bm{0)pB} = ^m(t) + «X™(i)/2. The general- 
ized time-dependent Redfield equation Q avoids the Markov 
approximation in the sense that the upper integration limit 
goes to t instead of oo. One also observes the usual oscil- 
lating terms that mix population and coherences. Next, we 
perform the secular approximation, essentially averaging over 
these fast oscillating terms. Here, we would like to study the 
effect of the Markovian versus NM decoherence rates and for- 
mulate the master equation such that the NMQJ method can be 
straightforwardly applied. The secular approximation is justi- 
fied in the slow decoherence regime when \lu — oj'\~^ td 
for all transition frequency differences, where t/j is a general 
time scale of decoherence. Finally, we assume that every chro- 
mophore is embedded in an identical phonon environment. 
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thus the TO subscript for the correlator and the response func- 
tion can be dropped |30|. One arrives at the master equation 
in the interaction picture: 



dt 



Ps 



(t) - ^*L(i,c.)[Aj„HA™H,p(i)] (7) 



+-f{t,u)A„,{uj)p{t)Al{uj) 



with the time-dependent Lamb shift. 



L{t,uj) = Im{ f dtie-'^^'^Citi)}, 
Jo 



and the time-dependent rates, 



7(i,w) = 2Re{ f dtie-"^*iC(ii)}. 
Jo 



(8) 



(9) 



A transformation into the Schrodinger picture can be readily 
performed, resulting in the usual system Hamiltonian commu- 
tator term of the master equation lfT4l : 



dt 



Psit) 



:[Hs + HLsit),p{t)] 



(10) 



muj 

Y,U{t,oj){Al{uj)A^{uj),p{t)}. 



The Hamiltonian HLs{t) = E™^ ^(^, 
leads to a Lamb-type renormalization of the energy levels. In 
the present work, we do not consider this term since we do not 
expect a qualitatively new contribution to the exciton dynam- 
ics 1,34. .36,1 . 



III. SPECTRAL DENSITY AND TIME-DEPENDENT 
RATES 

The main result for a Redfield master equation with- 
out Markov approximation is the time dependence of the 
rates. The decoherence rates depend on the phonon coupling 
strengths A^. The spectral density (units of frequency) de- 
scribes the coupling strength at a particular frequency: 



J{lo) 



(11) 



Assuming a continuous distribution of modes the spectral den- 
sity can be modeled with various functional forms. In molecu- 
lar energy transfer often an Ohmic spectral density with expo- 
nential or Drude-Lorentz cutoff is employed |24, 30, 38 39 1. 
In Ref. lfT3l non-Markovian phonon sidebands in fluorescence 
spectra of the B777 complex were reproduced with a super- 
Ohmic spectral density. In this paper, we assume an Ohmic 
spectral density with exponential cutoff: 



J(uj) — - — wexp 



(12) 



The relevant quantities are the cutoff ujc and the reorganiza- 
tion energy A. The cutoff determines the position of the peak 
of the spectral density and the reorganization energy is given 
hy X — h J duj^^^^. In Fig.jl (upper panel) we show the spec- 
tral density for a particular choice of parameters. We choose 
A = 30cm~^, which is typical for chromophores in photosyn- 
thetic systems |20|, and ojc ~ 30cm^^ which corresponds to 
relatively slow phonon modes. Note that typical transition fre- 
quencies in the single exciton manifold such as a; « 200cm^^ 
are located at the tail of the spectral density. The resulting 
Markovian relaxation rates are small. The strongly coupled, 
"off-resonant" modes at around 30cm^^ can lead to consider- 
able non-Markovian effects of the decoherence rates. 

For any spectral density and for the bosonic bath, the time- 
dependent decoherence rate is derived from Eq. (j9]): 

j{uj,t) = 2 / dCbJiCb) [n{Cd) — (13) 

Jo V uj + uj 



+{n{^) + 1) 



sin((cj — £l')t) 



Here, n{uj) is the bosonic distribution function. In the Marko- 
vian case {t oo), the spectral density is sampled only at 
the frequency us, seen from the limiting behavior of the terms 
^^j^ sin((ti; ± Oj)t). For the dephasing rate one obtains from 
Eq. ( pjj ) in the hmit a; ^ 0: 



Hlu \ sin(a)f) 



(14) 



In the Markovian limit, the dephasing rate becomes linearly 
proportional to the temperature and the derivative of the spec- 
tral density at zero frequency 1 14|. In the non-Markovian case, 
a greater part of the spectral density is taken into account. This 
can lead to rich behavior of both relaxation and dephasing 
rates. In Fig. [T] (lower panel), we show the rates that follow 
from the spectral density ( 12 1 and the above choice of param- 
eters. The relaxation rates in the NM case oscillate around 
the Markovian rates, have positive and negative regions, and 
finally converge to the Markovian rate on a time scale of Ips. 
The NM dephasing rate converges from below to the Marko- 
vian limit on a similar time scale. This transient regime has 
been discussed in terms of slippage in the initial conditions 
e.g. in L40il . 



IV. NON-MARKOVIAN QUANTUM JUMPS 

In this work, we perform a stochastic unraveling of the mas- 
ter equation with the non-Markovian quantum jump approach 
established in Ref. 1261 l27l . The master equation ([TOjl is pre- 
cisely in the form required for the NMQJ method. We give 
a brief summary of this technique. For every time t, one 
can separate the set of jump generators into A^{uj) 

and A^^ (uj) depending on the overall sign of the correspond- 
ing rate. That is for all yl+(a;) the rate is 'y^{uj,t) > 0, 
while for all A^^{llj) the rate is j^{u!,t) < 0. In the pres- 
ence of only positive channels, the original MCWF method 
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FIG. 1: Spectral density and resulting time-dependent decoherence 
rates from Eqs. l |13^ and l |14[ l. The upper panel shows the Ohmic 
spectral density with exponential cutoff for the parameters A = 
30cm~^ and ujc = 30cm^^. In the physical situation studied in the 
present work, the main strongly coupled modes are "off-resonant" to 
a transition frequency (200cm^^ in this figure). This leads to rich 
behavior of the corresponding rates at relevant time scales of around 
Ips. The relaxation rates oscillate, turn negative, and converge to 
their Markovian limit (lower panel, left); the inset shows a magnifi- 
cation at times just before Ips. The dephasing rate rises from zero 
and converges to the Markovian limit (lower panel, right). 



can be employed ||25l |26| . The particular structure of the 
jump generators (see previous sections) leads to a rel- 

atively straightforward description of the quantum mechanical 
ensemble. The density matrix at all time can be written as: 



pit) 



N 



\Ait)){Ait)\ 



^^|M)(M|. (15) 



M 



N 



Here, \ipi{t)) is the initial state with statistical weight 
iVo(t)/iV. The exciton states \AI) are as defined above and 
have a statistical weight NM{t)/N. Initially, iVo(O) = and 
at all times the numbers iVo(t) and Nmit) conserve proba- 
bility, i.e.No{t)+ J2m ^M{t) — N. The time evolution con- 
sists of propagation of \ipi{t)) and stochastic changes of the 
weights NQ{t) and NM{t)- In general, one defines the effec- 
tive Hamiltonian: 



ih 



Hes = Hs- - ^7(i,'^)4„H^m(^), (16) 



where the sum is over positive and negative channels. The 
NMQJ method now describes the time evolution of the ensem- 
ble p{t) as a wavefunction evolution of the ensemble states 



with Heft interrupted by probabilistic, discontinuous jumps 
corresponding to the jump operators of all channels. Consider 
now a particular ensemble member \ip{t)) at time t evolving 
for a small time step 6t. As in the MCWF, the no-jump evolu- 
tion is: 



1 



iSt 1 1 



\m) 



ii(i-fi7eff) mm' 



(17) 



The positive jumps occur with probability PmLo{t) — 
St j+{uj,t) (i}{t)\A+^uj)A+{Lu)\i}{t)) and an ensemble 
member jumps according to: 



(18) 



The negative jumps occur from the source state \ ip{t)) to target 
states \ip'{t)) if the source \ijj{t)) has the property that it can 
be reached by a jump with A" (w) from the target state: 



m)) = 



lll^™M^'W>ll 



\ij'{t + 6t)) (19) 



Note that a negative jump can "undo" positive jumps 
that occurred earlier. The negative jump probability de- 
pends on the target state \ip'{t)) and is Pmu){t) ~ 
^<5t|7+(c^,<)| (V'W|A^(^M™HI^'(i)>, where 7V'(t) 
is the number of ensemble members in the target state and 
N{t) are the number of ensemble members in the source state. 
A Monte Carlo unraveling according to this prescription is 
shown to be equivalent to master Eq. 



10 1, see Ref. |26|. 



Non-Markovian quantum jumps can explicitly lead to restored 
quantum coherence with this jump description that correctly 
handles negative rates in the master equation. 

We end this section with a note regarding positivity of the 



density matrix. The master equation ( 10 1 with time-dependent 
rates is not guaranteed to ensure positivity of the density ma- 
trix. However, the NMQJ method yields a simple criterion for 
detecting when positivity is about to be violated based on the 
negative jump probability |28 |. The P^nuji^) is inversely pro- 
portional to the number of ensemble members in the source 
state N{t). The case when N{t) becomes zero and the rate is 
negative at the same time is precisely when the master equa- 
tion would violate positivity. The interpretation of this is that 
the environment tries to "undo" an event that never happened. 
Thus, based on the singularity of the negative jump probabil- 
ity one can easily detect unphysical time evolution in the al- 
gorithm. All results in this work originate from physical time 
evolution. 



V. POPULATION BEATINGS IN A DIMER SYSTEM 

In this section, we show that oscillatory non-Markovian de- 
coherence rates can lead to beatings of site populations that 
do not occur in a Markovian treatment of the dynamics. The 
beatings arise from the coupling to slow modes in the envi- 
ronment. We discuss a dimer system consisting of two inter- 
acting chromophores in a structured phonon bath. The system 
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FIG. 2: Decoherence rates and density matrix elements for a dimer tiiat resembles site 1 and 2 of tlie Fenna-Mattiiews-Olson complex for 
Markovian and non-Markovian evolution. The left panels show dephasing (respective upper panel) and relaxation rate (respective lower 
panel) as a function of time. The time-dependent rates (blue) converge to their respective Markovian limits (red). Time evolution of the 
population elements (p^^ blue, P22 green) and coherence elements (Repj^j red, Impjj cyan) is displayed in the Markovian (middle panels) and 
non-Markovian (right panels) case for various temperatures. At room temperature NM population beatings are observed. 



Hamiltonian in the single exciton manifold is: 

Hs = e2|2)(2| + V12 (|1>(2| + |2)(1|) . (20) 

The eigenenergies are E2.1 — 62/2 ± 1 /2y^e^~+4y^ and the 

transition frequency is huj2i — {E2 ~ Ei) = -^/gj + 4V^2- 
The respective eigenstates are \Ei) — ci(l)|l) + C2(l)|2) 
and \E2) = ci(2)|l) + C2(2)|2), with ci(l) = -C2(2) = 
sm6 and ci(2) = 02(1) = cos 9 with the mixing angle 



tan 26* = 2Vi2/e2- The jump generators for relaxation are 
A2{uj) = -^i(t^) = l/2sin26'|£'i)(i;2| and their trans- 
pose conjugates. The jump generators for dephasing are 
^i(O) = sm^e\Ei){Ei\ + cos'^9\E2){E2\ and ^2(0) = 
cos^ 0\Ei){Ei\ + sin^ e\E2){E2\. 

For the simulations we take the excitonic Hamiltonian to 
be of a particular form: V12 = 87cm^^ and 62 = 120cm^^. 
This form is equal to the Hamiltonian of the site 1 and 2 sub- 
system in the Fenna-Matthews-Olson complex given in |24|. 
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Nevertheless, the effects presented here hold for a large vari- 
ety of dimer Hamiltonians. The initial state is localized at site 
1, i.e. p{0) ~ In Fig. [2] we show the time dependence 

of the decoherence rates (left panels) and the time evolution of 
the population and coherence elements of the density matrix 
in the site basis, i.e. PmnW = {m\p{t)\n). We compare dif- 
ferent temperatures in the Markovian (middle panels) and the 
non-Markovian description (right panels). The spectral den- 
sity parameters are reorganization energy A = 50cm^^ and 
cutoff u!c = 50cm^^. The rates (left panels in Fig. [2]) are sim- 
ilar as discussed for Fig. [T] The relaxation rates are generally 
smaller than the pure dephasing rate. The Markovian limit is 
reached on a time scale on the order of Ips. Higher temper- 
ature leads to larger oscillations of the decoherence rates. In 
the Markovian case (middle panels in Fig.|2]i, we observe that 
the population and coherence oscillations die out very fast, es- 
pecially at high temperatures. This is explained by the linear 
dependence of the dephasing rate on temperature. 

In the non-Markovian case, the dynamics is considerably 
different. At low temperatures the beatings are similar to the 
Markovian case; decoherence rates are generally rather small 
and differences in NM versus Markovian are not pronounced. 
At increasing temperature, the quantum mechanical beatings 
live slightly longer in the NM case due to the smaller dephas- 
ing rate at short times. Quantum beatings can be recognized 
by an oscillating imaginary part of the coherence element of 
the density matrix. At large temperatures, another type of 
beating arises, which is due to the oscillatory relaxation rates. 
It leads to beatings of the populations matrix elements and the 
real part of the coherence matrix element. These beatings can 
be interpreted as recurrence of the NM environment; energy is 
emitted from the system into the environment during the pos- 
itive regions of the decoherence rates and reabsorbed in the 
same decoherence channel during the negative regions. 



VI. TRANSPORT IN THE NON-MARKOVIAN REGIME 

In this section, we focus on transport properties in the non- 
Markovian regime. The general behavior of the exciton trans- 
port as a function of different parameters, contributions of 
various physical processes, and robustness of a chromophoric 
network can be investigated by theoretical measures such as 
the energy transfer efficiency and the transfer time ['36', '4T1. 
Trapping sites model the reaction centers where charge sepa- 
ration and energy storage occurs in the photosynthetic system, 
neglecting further chemical detail. In this paper, we utilize 
a simpler measure to elucidate energy transport. We define 
the integrated probability of a particular excitonic state up to 
a certain time r, the only free parameter. An explicit intro- 
duction of trapping sites and additional free parameters is not 
required. Formally, the measure is given by: 



Pi 



1 



M 



dt {M\pit)\M), 



(21) 



choose an exciton state \M) and focus on relaxation dynam- 
ics in the exciton basis. For example, in the Fenna-Matthews- 
Olson complex the exciton with the lowest energy, localized 
at site 3 and 4, would be an appropriate choice. A site-basis 
definition is straightforward within the NMQJ method. 



We analyze the transport properties of master equation (10 1 



and the NMQJ unraveling and compare to the Markovian 



regime. We choose the Hamiltonian parameters in Eq. ( 20 1 as 
V12 — 50cm^^ and 62 — 2Vi2- We assume that the system is 
initially in the energetically higher eigenstate | E2) and investi- 
gate relaxation to the lower lying eigenstate \Ei). We quantify 
the transport by the integrated probability of Eq. (21 1 using 
\Ei), i.e. Pi. In Fig.|3] we show the dependence of the trans- 
port on the essential parameters of the spectral density (reor- 
ganization energy A and cutoff coc) and the temperature. If not 
explored as variables, the default parameters are A — 30cm^^, 
u!c — 30cm~^, and room temperature (T = 300K). We inves- 
tigate two integration times, t = Ips and r = 4ps. The results 
are more pronounced for the shorter time. Shorter time scales 
are more relevant for smaller photosynthetic complexes such 
as the Fenna-Matthews-Olson complex |24l. As a result, we 
observe that transport can be enhanced in the NM situation 
where the rates are gives by Eq. ( 13 1. This study can be seen 



where r is the total integration time and \M) is a particular ex- 
citon state given by the problem at hand. For this measure we 



as an extension of ENAQT concepts to the non-Markovian 
regime. 

In the upper panels of Fig. [3] the dependence of Pi as a 
function of the reorganization energy A is investigated. Relax- 
ation and dephasing rates are scaled linearly with A. In gen- 
eral, the probability Pi increases as a function of the reorgani- 
zation energy. When relaxation rates are larger, thermal equi- 
libration of the exciton populations is faster At A = 30cm~^ 
and r = Ips the non-Markovian probability is Pi = 0.44 
while the Markovian probability is Pi — 0.27, a substantial 
difference for this rather small system. The improvement can 
be rationalized by the fact that the initially large NM relax- 
ation rates lead to fast equilibration. The negative regions 
of the rate partially "undo" the positive region but the inte- 
grated population of the target exciton is overall larger than 
in the Markovian case. For reorganization energies beyond 
«90cm^^ we observe positivity violating time evolution iden- 
tified with the criterion discussed earlier [281 . 

The middle panels of Fig.[3]show the dependence of Pi on 
the temperature. In the present case of energy transfer from 
a high to low exciton state, temperature can have an assisting 
effect for short times and for both Markovian and NM treat- 
ment |30|. For example, see the graphs for r = Ips. Increased 
thermal population of the phonon modes can lead to increased 
stimulated emission of exciton energy into the phonon bath 
and thus transport towards the lower exciton state. This effect 
becomes weaker for longer times, see the graphs for t = 4ps. 
Absorption of energy from the phonon bath comes into play 
which transports the excitation from the lower exciton state 
back to the higher one. The temperature is more significant 
in the NM regime since the Bose functions in the rate integral 
are sampled at all frequencies instead of only at uj2i- 

In the lower panels of Fig. [3] the averaged probability Pi 
is shown as a function of the cutoff utc of the spectral density. 
For both NM and Markovian cases, a larger cutoff and thus a 
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FIG. 3: The averaged probability measure Pi as defined in Eq. \2l\ as a function of the main decoherence parameters for a dimer system. 
Markovian (green) and non-Marlcovian (blue) cases are depicted for different temperatures and the central parameters of the spectral density 
(reorganization energy A, cutoff Uc)- Two integration times for the measure, r = Ips (left panels) and r = 4ps (right panels), are shown, the 
effects being more pronounced for the shorter time scale. The standard parameters are lUc = 30cm^^, A — 30cm^^, and room temperature 
T = 300K. 



stronger coupling of modes that are resonant with the transi- 
tion frequency leads to increased transport. The differences in 
terms of transport of both cases vanish when cUc ~ 60cm^^; 
resonant modes dominate the relaxation dynamics. However, 
in the presence of only slow modes, the NM treatment shows 
substantially larger transfer probabilities. For r = Ips and 
ujc — 20cm^^, we obtain Pi = 0.31 in the NM case and 
Pi = 0.06 in the Markovian case. This improvement is due 
to sampling of the a broader range of the spectral density in 
Eq. ( [T3] l. The physical interpretation is that the NM treat- 
ment allows the system to temporarily access energy non- 
conserving phonons for quantum jumps that would be inac- 
cessible otherwise. 



VII. FENNA-MATTHEWS-OLSON COMPLEX 

We also performed simulations for the Fenna-Matthews- 
Olson (FMO) complex. The FMO complex acts as an energy 
transfer wire in green sulphur bacteria Chlorobium tepidum 
(12 • It is subject of recent experimental ||23| and theoretical 
studies El [Ml [35] im. We derived the TCL master equa- 



tion ( 10 1 for the seven-site FMO subunit and performed sim- 
ulations with NMQJ method. We used the Hamiltonian of 
ref. | |24l and the spectral density ( 12 1 with A = 35cm^^ and 
LOc — 150cm^^ [30|. The decoherence rates for the 42 re- 
laxation channels (absorption + emission) and the 7 dephas- 
ing channels are evaluated. The rates oscillate and some have 
negative regions. The initial states are chosen to be localized 
at site 1 or 6, the sites that are close to the chlorophyll an- 
tenna complex. As a result, we find that the dynamics is not 
substantially affected by the time-dependent rates, see Fig.|4] 
Quantum beatings are slightly longer-lived for both tempera- 
tures 77K and BOOK and both initial states. This is because 
the NM dephasing rates converge from below to the Marko- 
vian limit similar to Fig. [T] (bottom right panel). The main 
relaxation rates stay positive and oscillate around their Marko- 
vian values. The spectral density is rather broad, covering all 
transition frequencies, cf. Fig. 2 of [35], such that the effects 
described in the previous sections turn out to be not dominant. 
The Markovian approximation alone in the presence of the 
other approximations such as Born and secular does not have a 
substantial effect. Recently, Ishizaki and Fleming utilized the 
hierarchical equation of motion approach for explaining long- 
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FIG. 4: Population of sites in the Fenna-Matthews-Olson complex for Markovian (dashed) and non-Markovian (solid) cases. Sites that are 
close to the chlorophyll antenna are taken to be the initial states, i.e. site 1 (upper panel) or site 6 (lower panel). Parameters are ujc ~ 150cm~^, 
A = 35cm" \ and T = 300K. 



lived coherences in the FMO complex ll24l . Since Born and 
secular approximations are avoided for Gaussian fluctuations, 
this approach has a larger range of validity than the Redfield 
model, especially at large temperatures, and correctly incor- 
porates molecular reorganization effects. 

Vm. CONCLUSION 

In conclusion, we have applied the non-Markovian quan- 
tum jump method to excitonic energy transfer. The NM de- 
coherence rates resulting from a time-convolutionless treat- 
ment of the master equation are oscillatory and negative for 
parameter regimes and time scales that are relevant to the 
problem. In the present work, NM effects are large when a 
system is strongly coupled to "off-resonant" phonon modes 
of the environment. These slow modes can lead to popula- 
tion beatings at room temperature, which are a signature of 
bath recurrence effects. Additionally, our computations show 
that Markovian versus non-Markovian dynamics can thus cru- 
cially affect transport dynamics. Quantum transport can be 
enhanced over strictly Markovian dynamics due to a sam- 
pling of broader regions of the spectral density. We thus have 
provided a non-Markovian extension to recent environment- 



assisted quantum transport (ENAQT) concepts. For example, 
a system with a transition of around 140cm" ^ shows consid- 
erable NM improvement of transport in the presence of strong 
modes at around 30cm~^. 

Recently, Jang et al. ifTTI developed a novel theory of co- 
herent resonance energy transfer. A small polaron transforma- 
tion is applied before the second-order time-convolutionless 
expansion that leads to time-dependent decoherence rates. 
Nonequilibrium reorganization effects are taken into account 
by the exciton-phonon dressed state description. This treat- 
ment can lead to an increased range of validity with respect 
to the exciton-phonon couplings and temperatures compared 
to the standard Redfield approach. In this context, the NMQJ 
method in its present form or with suitable extensions could 
prove especially powerful to efficiently simulate larger donor- 
acceptor systems and to correctly incorporate negative deco- 
herence rates in a quantum jump description. Fuithermore, 
the NMQJ method can be used for the stochastic computation 
of spectroscopic signals Il43 1 i44l . 
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